;;------------------------------------------------------------------------- 
;;   - Drawing box plots
;;   - Using a box plot to show the the median, minimum/maximum value, and the 25th/75th percentiles of two timeseries
;;   - Adding text to a box plot
;;   - Sorting data
;;   - Setting the mininum/maximum value of the Y axis in a box plot
;;   - Generating dummy data using "random_normal"
;;------------------------------------------------------------------------- 
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
load "$KAILIBROOT/lib_radon.ncl"
load "$KAILIBROOT/lib_boxplot.ncl"  ;; use my own boxplot 

begin

  plottype = "pdf"
  plottype = "eps"
  plotname = "fig_box_" 


  USA087   = False
  ionprod  = False  

  PP = "02"  
  PP = "01"  
  PP = "04"  
  PP = "03"  

  load "load_PP.ncl" 

  StringFontHeightF  = 0.035 
  cnLineLabelFontHeightF = 0.025 
  tmXBLabelFontHeightF = 0.022
  tmYLLabelFontHeightF = 0.022
  tiMainFontHeightF  = 0.022
  tiXAxisFontHeightF = 0.022
  tiYAxisFontHeightF = 0.022
  lbLabelFontHeightF = 0.030 

  cpp = "" 

  sname = asciiread("radon_sitename.dat",-1,"string")

  ns = dimsizes(sname)
  np = 24

  plot = new(ns,graphic)  
  plop = new(np,graphic)  


;;------------------------------------------------------------------------- 
;; plot index 
;;------------------------------------------------------------------------- 
  idp = new((/np/),"integer") 

  idp(0) = ind(sname.eq."freiburg")
  idp(1) = ind(sname.eq."schauinsland")
  idp(2) = ind(sname.eq."hohenpeissenberg")
  idp(3) = ind(sname.eq."gif_sur_yvette")
  idp(4) = ind(sname.eq."zingst")
  idp(5) = ind(sname.eq."lutjewad")

  idp(6) = ind(sname.eq."hongkong")
  idp(7) = ind(sname.eq."beijing")
  idp(8) = ind(sname.eq."wuhan")
  idp(9) = ind(sname.eq."changchun")
  idp(10)= ind(sname.eq."xi_an")
  idp(11)= ind(sname.eq."huhehaote") 

  idp(12)= ind(sname.eq."cincinnati")
  idp(13)= ind(sname.eq."bermuda")
  idp(14)= ind(sname.eq."socorro")
 ;idp(14)= ind(sname.eq."grifton0001")
  idp(15)= ind(sname.eq."para")
  idp(16)= ind(sname.eq."cape_grim")
  idp(17)= ind(sname.eq."cape_point") 

  idp(18)= ind(sname.eq."amsterdam_island")
  idp(19)= ind(sname.eq."crozet")
  idp(20)= ind(sname.eq."kerguelen")
 ;idp(21)= ind(sname.eq."macquarie_island")
  idp(21)= ind(sname.eq."mauna_loa")
  idp(22)= ind(sname.eq."dumont")
  idp(23)= ind(sname.eq."mawson") 


  str = new(np,"string") 
 
  str(0 ) = "Freiburg" 
  str(1 ) = "Schauinsland" 
  str(2 ) = "Hohenpeissenberg" 
  str(3 ) = "Gif-sur-Yvette" 
  str(4 ) = "Zingst" 
  str(5 ) = "Lutjewad" 

  str(6 ) = "Hong Kong" 
  str(7 ) = "Beijing" 
  str(8 ) = "Wuhan" 
  str(9 ) = "Changchun" 
  str(10) = "Xi'an" 
  str(11) = "Huhehaote"  

  str(12) = "Cincinnati" 
  str(13) = "Bermuda" 
  str(14) = "Socorro" 
 ;str(14) = "Grifton" 
  str(15) = "Para" 
  str(16) = "Cape Grim" 
  str(17) = "Cape Point" 

  str(18) = "Amsterdam Island" 
  str(19) = "Crozet" 
  str(20) = "Kerguelen" 
 ;str(21) = "Macquarie Island" 
  str(21) = "Mauna Loa" 
  str(22) = "Dumont" 
  str(23) = "Mawson" 


  print("idp: "+idp)  

  obsfilepath = "/pf/m/m222044/post/my_script/radon_cons_monthly/dir.obs/data/"  

 
;;------------------------------------------------------------------------- 
;; load settings by user 
;;------------------------------------------------------------------------- 


;;------------------------------------------------------------------------- 
;; filename 
;;------------------------------------------------------------------------- 
  filenama = "box_data_ANN.nc" 
  filenamb = "box_data_DJF.nc" 
  filenamc = "box_data_MAM.nc"
  filenamd = "box_data_JJA.nc"
  filename = "box_data_SON.nc"

  haveobs = isfilepresent(filenama) 

  fla = addfile(filenama,"r") 
  flb = addfile(filenamb,"r") 
  flc = addfile(filenamc,"r") 
  fld = addfile(filenamd,"r") 
  fle = addfile(filename,"r") 

;;------------------------------------------------------------------------- 
;; read data 
;;------------------------------------------------------------------------- 
  ppa = fla->boxa ;; percentile 
  ppb = flb->boxa
  ppc = flc->boxa
  ppd = fld->boxa
  ppe = fle->boxa

  paa = fla->boxb ;; mean 
  pab = flb->boxb
  pac = flc->boxb
  pad = fld->boxb
  pae = fle->boxb

  pba = fla->boxc ;; mean 
  pbb = flb->boxc
  pbc = flc->boxc
  pbd = fld->boxc
  pbe = fle->boxc

  is = 0 

  print("") 
  print("read data") 
  print("") 


;;------------------------------------------------------------------------- 
;; merge data  
;;------------------------------------------------------------------------- 

  nep = 6 ;; number of exp, available model data 
  npp = 4 ;; number of exp, plot 
  npt = 5 
  nsn = 5 

  xpp = new((/ns,nep,nsn,npt/),"float") 
  xpa = new((/ns,nep,nsn/),"float") 
  xpb = new((/ns,nsn/),"float") 

  nsn = 5 
 
  do isn = 0, nsn-1
     xpp(:,:,0,:) = ppa(:,:,:)  
     xpp(:,:,1,:) = ppb(:,:,:)  
     xpp(:,:,2,:) = ppc(:,:,:)  
     xpp(:,:,3,:) = ppd(:,:,:)  
     xpp(:,:,4,:) = ppe(:,:,:)  
  end do 
 
  xpa(:,:,0) = paa(:,:)  
  xpa(:,:,1) = pab(:,:)  
  xpa(:,:,2) = pac(:,:)  
  xpa(:,:,3) = pad(:,:)  
  xpa(:,:,4) = pae(:,:)  

  xpb(:,0) = pba(:)  
  xpb(:,1) = pbb(:)  
  xpb(:,2) = pbc(:)  
  xpb(:,3) = pbd(:)  
  xpb(:,4) = pbe(:)  

;;------------------------------------------------------------------------- 
;; update plot data 
;;------------------------------------------------------------------------- 
  load "plotbox_panel_boxval.ncl" 

;;------------------------------------------------------------------------- 
;; reshape array 
;;------------------------------------------------------------------------- 

 ;pp = new((/ns,nep*nsn,npt/),"float") 
 ;pa = new((/ns,nep*nsn/),"float") 
  pp = new((/ns,npp*nsn,npt/),"float") 
  pa = new((/ns,npp*nsn/),"float") 

  ;printVarSummary(pp) 
  ;printVarSummary(xpp)  

  isp = 0 

  do isn = 0, nsn-1 ;; season second 
  do iep = 0, npp-1   ;; exp first  
     ;;print("isn,iep,isp: "+isn+" "+iep+" "+isp) 
     ;;       isn,iep,isp: 4 3 20
     pp(:,isp,:) = xpp(:,iep,isn,:) 
     pa(:,isp)   = xpa(:,iep,isn) 
     isp = isp + 1  
  end do 
  end do 

;;------------------------------------------------------------------------- 
;; radon activity to ion-production rate 
;;------------------------------------------------------------------------- 

  if(ionprod) then 
     pp = RN222_ionprod(pp)  
     pa = RN222_ionprod(pa)  
  end if  
 
  pp = where(pp.lt.1.E-06,-10000,pp) 
  pa = where(pa.lt.1.E-06,-10000,pa)   


;;------------------------------------------------------------------------- 
;; for clim show 
;;------------------------------------------------------------------------- 
  clim = True

  if (clim) then
     xpp(:,2,:,:) = xpp(:,3,:,:)
     xpa(:,2,:)   = xpa(:,3,:)
     xpp(:,3,:,:) = xpp(:,5,:,:)
     xpa(:,3,:)   = xpa(:,5,:)
  end if

;;------------------------------------------------------------------------- 
;; define x-axis 
;;-------------------------------------------------------------------------  
  ppx = new( npp*nsn, float )

  do isn=0,nsn-1
  do iep=0,npp-1
      ppx(isn*npp+iep) = 5 + isn*(nsn*0.6) + iep*0.4
  end do
  end do 

  print("") 
  print("compute vars") 
  print("") 





;;------------------------------------------------------------------------- 
;; open wks 
;;------------------------------------------------------------------------- 

  wks = gsn_open_wks(plottype,plotname+PP) 

  res = True 

;;------------------------------------------------------------------------- 
;; change plot size 
;;------------------------------------------------------------------------- 
  res@mpShapeMode  = "FreeAspect"
  res@vpWidthF      = 0.7
  res@vpHeightF     = 0.35


;;------------------------------------------------------------------------- 
;; strings 
;;------------------------------------------------------------------------- 

 ;res@tmXBLabels = (/"OBS","E0","E1","E2","E3"/) ; labels for each box
  res@tmXBLabels = (/"OBS","E0","E2","E3", \
                     "OBS","E0","E2","E3", \ 
                     "OBS","E0","E2","E3", \
                     "OBS","E0","E2","E3", \
                     "OBS","E0","E2","E3"/) ; labels for each box
  res@tmXBLabels = (/"","    ALL","","", \
                     "","    DJF","","", \
                     "","    MAM","","", \
                     "","    JJA","","", \
                     "","    SON","",""  /)

 ;res@tiMainString     = sname(is) ;;+ ", "+season 
 ;res@tiXAxisString    = ""

  if(ionprod) then 
  res@tiYAxisString    = " ~C~ IPRR (cm~S~-3~N~s~S~-1~N~) "            ; y-axis label    
  else
  res@tiYAxisString    = "Radon activity (mBq m~S~-3~N~STP) "    ; y-axis label  
  res@tiYAxisString    = "Radon (mBq m~S~-3~N~STP) "    ; y-axis label     
  end if 
;;----------------------------------------------------------------------
;; tickmark 
;;----------------------------------------------------------------------
;;res@tmYLMode = "Explicit"	
;;res@tmYLValues = (/100.,200.,400.,600.,800.,1000.,2000.,3000.,4000./)
;;res@tmYLLabels = (/"100","200","400","600","800","1000","2000","3000","4000"/)
;;res@tmYLMinorValues  = ispan(0,4000,200)

  ;res@tmYROn = True 
  ;res@tmXTOn = True 
  ;res@tmYRBorderOn = True 
  ;res@tmXTBorderOn = True 



  res@tmXBLabelFontHeightF = tmXBLabelFontHeightF
  res@tmYLLabelFontHeightF = tmYLLabelFontHeightF

;;----------------------------------------------------------------------
;; axis  
;;----------------------------------------------------------------------
  res@tiXAxisFontHeightF = tiXAxisFontHeightF
  res@tiYAxisFontHeightF = tiYAxisFontHeightF

;;------------------------------------------------------------------------- 
;; legend 
;;------------------------------------------------------------------------- 
 ;res@pmLegendDisplayMode    = "Always"              ; turn on legend
 ;res@pmLegendSide           = "Top"                 ; Change location of 
 ;res@pmLegendParallelPosF   = .90                   ; move units right
 ;res@pmLegendOrthogonalPosF = -0.8                  ; more neg = down
 ;res@pmLegendWidthF         = 0.12                  ; Change width and
 ;res@pmLegendHeightF        = 0.25                  ; height of legend.
 ;res@lgLabelFontHeightF     = .02                   ; change font height
 ;res@lgPerimOn              = False                 ; no box around
 ;res@xyExplicitLegendLabels = (/"20N","30N","40N","50N"/)


;;------------------------------------------------------------------------- 
;; resources for polylines that draws the boxes
;;------------------------------------------------------------------------- 
  resl                   = True                        
  resl@gsLineThicknessF  = 1.0                 ; line thickness 

;;------------------------------------------------------------------------- 
;; resources that control color and width of boxes
;;------------------------------------------------------------------------- 
  resb          = True                  
  resb@boxWidth = .30                           ; Width of box (x units)
  resb@boxColors = (/"black","blue","red","green", \  
                     "black","blue","red","green", \
                     "black","blue","red","green", \ 
                     "black","blue","red","green", \
                     "black","blue","red","green"  /)     ; Color of box(es)

;;------------------------------------------------------------------------- 
;; tickmark 
;;------------------------------------------------------------------------- 
  res@tmXBMajorOutwardLengthF = 0. ;;0.021   ; set tickmarks 
  res@tmXBMinorOutwardLengthF = 0. ;;0.015
  res@tmXBMajorLengthF = 0.006
  res@tmXBMinorLengthF = 0.006
  res@tmYLMajorLengthF = 0.010
  res@tmYLMinorLengthF = 0.006


  print("") 
  print("making plot") 
  print("") 


;;------------------------------------------------------------------------- 
;; loop over site
;;------------------------------------------------------------------------- 

  do is = 0,ns-1 

;;------------------------------------------------------------------------- 
;; only loop over wanted index 
;;------------------------------------------------------------------------- 

  if (any(is.eq.idp)) then 

;;------------------------------------------------------------------------- 
;; min & max values for y-axis 
;;------------------------------------------------------------------------- 

  print("set min & max values") 

  res@trXMinF = 3.5  
  res@trXMaxF = 20

  res@trYMinF = 0.

  ismax = max(pp(is,:,:)) .gt. max(pa(is,:)) 

  if(ismax) then  
     res@trYMaxF = max(pp(is,:,:)) * 1.2 
  else
     res@trYMaxF = max(pa(is,:))   * 1.2 
  end if  


;;------------------------------------------------------------------------- 
;; set max & min value 
;;------------------------------------------------------------------------- 
  load "plotbox_panel_maxmin.ncl" 

  ;res@tmYRMode = "Explicit"
  ;res@tmYRLabelsOn = True  
  ;res@tmYRLabels = RN222_ionprod_conc(tmYLValues)  


;;------------------------------------------------------------------------- 
;; do not print axis title 
;;------------------------------------------------------------------------- 
  ;if (is.eq.idp(1)) then 
  ;   res@tiYAxisString = ""  
  ;end if 

  ;if (is.eq.idp(3)) then 
  ;   res@tiYAxisString = ""  
  ;end if 

  ;if (is.eq.idp(5)) then 
  ;   res@tiYAxisString = ""  
  ;end if 

;;------------------------------------------------------------------------- 
;; make plot 
;;------------------------------------------------------------------------- 

  do ip = 0,np-1
  if(is.eq.idp(ip)) then 
     res@tiMainString = str(ip)
     res@tiMainFontHeightF = tiMainFontHeightF
  end if 
  end do 


  plot(is) = boxplot2(wks,ppx,pp(is,:,:),resb,res,resl)      
 ;plot(is) = boxplot3(wks,ppx,pp(is,:,:),resb,res,resl)      

;;------------------------------------------------------------------------- 
;; adding additional text to the plot
;;------------------------------------------------------------------------- 
;  tres = True
;  str  = "E0: WCRP1995 ~C~"+\
;         "E2: Schery (1998) scaled by 1.6 ~C~"+ \
;         "E3: Merged radon flux map" 
;         ;Szegvary et al. (2007)
;         ;Zhuo et al. (2008)
;         ;Griffiths et al. (2010)
;         ;Schery (1998)
;  tres@txFontHeightF = 0.015
;  gsn_text_ndc(wks,str,0.60,0.73,tres)



  end if ;; if plot index 

;;------------------------------------------------------------------------- 
;; loop over site end 
;;------------------------------------------------------------------------- 

  end do ;; ns 




  res@gsnDraw  = False                          ; don't draw
  res@gsnFrame = False                          ; don't advance frame

;;------------------------------------------------------------------------- 
;; add some polymarkers
;;------------------------------------------------------------------------- 
  print("") 
  print("making markers") 
  print("") 

  mres               = True                     ; marker mods desired
  mres@gsMarkerIndex = 1                        ; polymarker style
  mres@gsMarkerSizeF = 20.                      ; polymarker size
  mres@gsMarkerColor = "black"                  ; polymarker color

  dum = new((/ns,nsn*nep/),graphic) 

  do ip = 0,np-1 
     plop(ip) = plot(idp(ip)) 
     do isp = 0,nsn*npp-1 
        if(isp.eq.0 .or. isp.eq.4 .or. isp.eq.8 .or. isp.eq.12 .or. isp.eq.16) then 
           mres@gsMarkerColor = "black" 
        end if 
        if(isp.eq.1 .or. isp.eq.5 .or. isp.eq.9 .or. isp.eq.13 .or. isp.eq.17) then 
           mres@gsMarkerColor = "blue" 
        end if 
        if(isp.eq.2 .or. isp.eq.6 .or. isp.eq.10.or. isp.eq.14 .or. isp.eq.18) then 
           mres@gsMarkerColor = "red" 
        end if 
        if(isp.eq.3 .or. isp.eq.7 .or. isp.eq.11.or. isp.eq.15 .or. isp.eq.19) then 
           mres@gsMarkerColor = "green" 
        end if 
        dum(idp(ip),isp) = gsn_add_polymarker(wks,plop(ip),ppx(isp),pa(idp(ip),isp),mres) 
     end do 
  end do 

;;------------------------------------------------------------------------- 
; adding additional text to the plot
;;------------------------------------------------------------------------- 
  tres               = True                    ; text mods desired 
  tres@txFontHeightF = 0.012

  if(False) then 
     load "plotbox_panel_ndc.ncl" 
     ;drawNDCGrid(wks) 
  end if ;; if(False) 

;;----------------------------------------------------------------------
;; total plot index 
;;----------------------------------------------------------------------

  if(PP.eq."01") then  
  ii = 0  
  end if 

  if(PP.eq."02") then  
  ii = 6  
  end if 
  
  if(PP.eq."03") then   
  ii = 12
  end if 

  if(PP.eq."04") then   
  ii = 18 
  end if 

;;----------------------------------------------------------------------
;; panel plot
;;----------------------------------------------------------------------

  resp                    = True   ; panel mods desired
  resp@gsnFrame           = False  ; don't advance frame yet
 ;resp@gsnPanelBottom     = 0.0    ; space for label bar
 ;resp@gsnPanelTop        = 0.5    ; only panel on lower half of page
 ;resp@gsnPanelLabelBar   = True   ; common label bar
 ;resp@pmLabelBarWidthF   = 0.8    ; label bar width

  resp@gsnPanelYWhiteSpacePercent = 5
  resp@gsnPanelXWhiteSpacePercent = 2
 
  gsn_panel(wks,plop(ii:ii+5),(/3,2/),resp)

  system("ps2pdf "+plotname+PP+".eps") 

end






